function [parameters_est]=estimate_parameters(cfg)
% function [parameters_est]=estimate_parameters(cfg)
% 
% 
% % Assuming you've created the system.txt file to 
% % describe the system, this makes sure all the 
% % files are up to date:
% build_system;
% 
% % here we define the parameters to estimate:
% % the order is important: first system then variance
% to_estimate = {'Vp'
%                'kel'}';
% % get the system information:
% cfg = auto_fetch_system_information;
% % define the parameter set and the parameters
% % we want to include in the estimation
% cfg = select_set(cfg, 'default', to_estimate);
% 
% % if you want to estimate all  the parameters then 
% % select the parameter set in the following way:
% % cfg = select_set(cfg, 'default')
% 
% % You have to create the observation_details.m file: 
% % to see a description of the expected output type:
% % help check_cfg
% cfg.estimation.observation_function = 'observation_details';
%
% % This will display the final estimates and
% % text to update the system.txt file
% % set it to 'no' if you don't want to see anything
% cfg.options.verbose                 = 'yes';
%
%
% % The 'effort' field tells the estimation routine to try "harder" to find a
% % good parameter estimate. By default it will be assumed to be 1, and will just 
% % estimate the parameters like normal. 
%
% cfg.estimation.effort = 1;
%
% % Increasing this number will implement a sort of simulated annealing where 
% % moving around the parameter space adding noise to the current estimates in 
% % an attempt to find a better solution. The larger this number, the longer
% % it will take.
%
% % see the help for optimset for the different options
% % but this is how they are included with a couple examples:
% cfg.estimation.options = optimset('Display',   'iter', ...
%                                   'MaxIter',    3000,  ...
%                                   'MaxFunEval', 3000);
%

  if( isfield(cfg.estimation, 'options'))
    options = cfg.estimation.options  ;
  else
    options = optimset();
  end

 
  if(~isfield(cfg.estimation, 'effort'))
    cfg.estimation.effort = 1;
  end

%
% If there are variance parameters make sure the order is correct
%

%
%
% JMH add checks in check_cfg for ml vs wls
%
% if wls see if there are any variance parameters
% selected for estimation
%
% if ml make sure there are variance parameters
%
% run check_cfg to make sure it's all setup
%
% JMH modify output to update system.txt in order to include the num2text
% stuff in generate_report.m
%
%


  run_estimation = 'yes';
  if(~strcmp(check_cfg(cfg), 'good'))
    disp('#-> Check your cfg data structure');
    run_estimation = 'no';
  end

  if(~isfinite( calculate_objective(cfg.estimation.parameters.guess, cfg)))
    disp('#-> The objective function evaluates as either Inf or NaN');
    disp('#-> for the initial guess. Check your initial parameter guess');
    disp('#-> or the observation_details function');
    run_estimation = 'no';
  end


  if(strcmp(run_estimation, 'yes'))
    % initializing the best objective function value to 
    % infinity, this will be updated below
    best_obj      = inf;
    best_params   =  cfg.estimation.parameters.guess ;
    initial_guess = best_params;
    vp(cfg,  '-----------------------------------');
    vp(cfg,  'Beginning parameter estiamtion');
    if(strcmp(cfg.estimation.objective_type, 'wls'))
      vp(cfg,  '   Obj: Weighted Least Squares');
    elseif(strcmp(cfg.estimation.objective_type, 'ml'))
      vp(cfg,  '   Obj: Maximum Likelihood');
    end

    if(cfg.estimation.effort >1)
      vp(cfg,  'Pseudo simulated annealing selected:');
      vp(cfg,  '   Sit back and have a cup of coffee;');
      vp(cfg,  '   this is going to take a while');
    end
    vp(cfg,  '-----------------------------------');

    for est_ctr = 1:cfg.estimation.effort
       if(cfg.estimation.effort >1)
         if(est_ctr/cfg.estimation.effort <= .2)
           options =  optimset('Display',   'none', ...
                               'maxfuneval', 300, ...
                               'maxiter',    300);
           sigma     = .5;
           stage_ctr = 1;
         elseif(est_ctr/cfg.estimation.effort <= .5)
           options =  optimset('Display',   'none', ...
                               'maxfuneval', 300, ...
                               'maxiter',    300);
           sigma     = .4;
           stage_ctr = 2;
         elseif(est_ctr/cfg.estimation.effort <= .7)
           options =  optimset('Display',   'none', ...
                               'maxfuneval', 200, ...
                               'maxiter',    200);
           sigma     = .20;
           stage_ctr = 3;
         elseif(est_ctr/cfg.estimation.effort < 1)
           options =  optimset('Display',   'none', ...
                               'maxfuneval', 200, ...
                               'maxiter',    200);
           sigma     = .10 ;
           stage_ctr = 3;
         else
           options =  optimset('Display',   'Iter', ...
                               'maxfuneval', 500, ...
                               'maxiter',    500);
           sigma     = .0001;
           stage_ctr = 4;
         end
         
         vp(cfg,sprintf('Stage %d (sigma = %.1e) iteration %d of %d (%s)', ...
                         stage_ctr,...
                         sigma,...
                         est_ctr,...
                         cfg.estimation.effort, ...
                         datestr(now)));
           
       else
           sigma     = .0001;
       end


      
       % not trying too hard, just running fminsearch:
       [current_est, current_obj ] = ...
       fminsearch('calculate_objective', ...
                  initial_guess,         ...
                  options,               ...
                  cfg);
       

       % disp(sprintf('Intermediate Obj: %.6e', current_obj));

       % rebounding the parameters
       current_est    = bound_parameters(current_est   , ...
                                         cfg.estimation.parameters.lower_bound, ...
                                         cfg.estimation.parameters.upper_bound);
       
       % if we found a solution better than the previous best,
       % then we update our best parameter set and objective value
       if(current_obj < best_obj)
         best_params = current_est;
         best_obj    = current_obj;
         disp(sprintf('New Objective found: %.6e', best_obj));
       end

       % Adding noise to the best estimate so far  but only if 
       % we're not at the last stage that one starts with the
       % best guess
       if(est_ctr < cfg.estimation.effort - 1)


         initial_guess  = add_noise(best_params,  sigma.*ones(size(initial_guess)), 'log_normal');
         % rebounding guess
         initial_guess  = bound_parameters(initial_guess , ...
                                           cfg.estimation.parameters.lower_bound, ...
                                           cfg.estimation.parameters.upper_bound);
       else
         % On the last stage we start with whatever our 
         % best solution is so far, so best_params is stored
         % in initial_guess
         vp(cfg, 'On the last stage');
         initial_guess = best_params;
       end

    end

    parameters_est = best_params;
   
   
    try
      % We're done with the estimation, now we try to calculate the solution
      % statistics. This can be numerically tricky, so we have a fail
      % condition below:
      save /tmp/ss.mat;
      calculate_objective(parameters_est, cfg)
      ss = solution_statistics(parameters_est, cfg);
      generate_report(parameters_est, ss, cfg);
      if(strcmp(cfg.options.verbose, 'yes'))
        type('report.txt');
      end
   
      if(strcmp(cfg.options.verbose, 'yes'))
        vp(cfg, 'If you''re happy with the results, the following');
        vp(cfg, 'can be used to update system.txt file. Just copy, ');
        vp(cfg, 'paste, and delete the previous entries');
        for parameter_idx =1:length(parameters_est)
          ptmp.set_name  = cfg.estimation.set_name;
          ptmp.value     = var2string(parameters_est(parameter_idx), 12);
          ptmp.name      = cfg.estimation.parameters.names{parameter_idx};
          ptmp.units     = cfg.estimation.parameters.units{parameter_idx};
          ptmp.grouping  = cfg.parameters.matrix{(strcmp(cfg.parameters.matrix(:,1) , cfg.estimation.parameters.names(parameter_idx))), end};
          ptmp.lb_number = cfg.estimation.parameters.lower_bound(parameter_idx);
          ptmp.editable  = cfg.estimation.parameters.editable{parameter_idx};
          ptmp.ub_number = cfg.estimation.parameters.upper_bound(parameter_idx);
        
          ptmp.ub_text   = num2str(ptmp.ub_number);
      
          if(strcmp(cfg.estimation.set_name, 'default'))
            if(ptmp.lb_number == eps)
              ptmp.lb_text   = 'eps';
            else
              ptmp.lb_text   = num2str(ptmp.lb_number);
            end
        
        
            disp(sprintf('<P> %-10s  %-10s %8s %10s   %-8s %5s %10s ', ...
                          ptmp.name,     ...
                          ptmp.value,    ...
                          ptmp.lb_text,  ...
                          ptmp.ub_text,  ...
                          ptmp.units,    ...
                          ptmp.editable, ...
                          ptmp.grouping));
          else
            disp(sprintf('<PSET:%s:%s> %-1s', ...
                         ptmp.set_name,   ...
                         ptmp.name, ...
                         ptmp.value));
          end
        end
      end
    catch
      % If we fail to calculate the solution statistics
      % a message is displayed and the parameter values are
      % also given back.
      disp('#--> Solution statistics calculation failed');
      disp('#--> This can happen when you have a parameter');
      disp('#--> set that makes the system stiff. Below are');
      disp('#--> the final parameter estimates:');
      for parameter_idx =1:length(parameters_est)
        disp(sprintf('#--> %s = %.4f', ...
               cfg.estimation.parameters.names{parameter_idx}, ...
               parameters_est(parameter_idx)));
      end
    end

  else
  % creating an emapty parameters_est 
  % to return 
  parameters_est = [];
  end


